suppressPackageStartupMessages(library(scater))
suppressPackageStartupMessages(library(Seurat))
data_path = "D:/Dropbox/매뉴얼/SC_tutorial"
figure_path = "D:/Dropbox/매뉴얼/SC_tutorial"
setwd(data_path)
Load scesetFiltered, hvg
scesetFiltered = readRDS( file= "scsetFiltered_norm.rds")
load(file="hvg.RData")
Create seurat object
seuset <- CreateSeuratObject(
raw.data = counts(scesetFiltered)
)
seuset@data = exprs(scesetFiltered)
seuset@scale.data = exprs(scesetFiltered)
seuset@var.genes = rownames(hvg)
PCA
PCA = 50
seuset <- RunPCA(seuset, pcs.compute = PCA, weight.by.var = FALSE, do.print = TRUE)
## [1] "PC1"
## [1] "Hbb-bt" "Opalin" "Vtn" "Mog" "Trem2" "Mag" "Selplg"
## [8] "Aif1" "Mobp" "Gsn" "Mal" "Fcrls" "Tmem119" "Fcer1g"
## [15] "Csf1r" "Ctla2a" "Rgs5" "Flt1" "Cldn11" "Hbb-bs" "Pglyrp1"
## [22] "Hmgb2" "Ptgds" "Tyrobp" "Fgfr3" "Ly6c1" "Fam107a" "Cldn10"
## [29] "P2ry12" "Id1"
## [1] ""
## [1] "Malat1" "Fth1" "Tmsb4x" "mt-Nd1" "Calm2" "mt-Cytb" "Ftl1"
## [8] "Sepw1" "Snca" "Rps5" "Rps9" "Rplp0" "Ubb" "Pcp4"
## [15] "Rps4x" "Ppp3ca" "Meg3" "Rpl10" "Ptma" "Aldoa" "Ybx1"
## [22] "Tuba1a" "Zbtb20" "Bex2" "Tmsb10" "mt-Nd2" "Mllt11" "Ypel3"
## [29] "Mt3" "Rps29"
## [1] ""
## [1] ""
## [1] "PC2"
## [1] "Cst3" "Apoe" "Dbi" "Cd81" "Ptn" "Aldoc"
## [7] "Glul" "Atp1a2" "Slc1a2" "Fabp7" "Mt1" "Clu"
## [13] "Cpe" "Ftl1" "Itm2b" "Sparcl1" "Pantr1" "Tsc22d4"
## [19] "Slc1a3" "Ndrg2" "Cspg5" "Ckb" "Lamp1" "Scd2"
## [25] "S100a16" "Car2" "Plpp3" "Prdx6" "Cd63" "Serpine2"
## [1] ""
## [1] "Snca" "Pcp4" "Meg3" "Ppp3ca" "Calm2" "Ncdn" "Olfm1"
## [8] "Camk2b" "Bex2" "Sncb" "Mllt11" "Cplx2" "Cnih2" "Snrpn"
## [15] "Ndrg4" "Ctxn1" "Chn1" "Map1b" "Snap25" "Ly6h" "Ywhah"
## [22] "Camk2a" "Ncald" "Plppr4" "Rtn1" "Atp1b1" "Grin2b" "Nsf"
## [29] "Caly" "Atp2b1"
## [1] ""
## [1] ""
## [1] "PC3"
## [1] "Mt3" "Mt1" "Slc1a2" "Tspan7" "Chchd10" "Ncdn" "Atp1b1"
## [8] "Sparcl1" "Aldoc" "Ppp3ca" "Aldoa" "Chn1" "Apoe" "Ppp1r1a"
## [15] "Gpm6a" "Clu" "Camk2a" "Atp1a2" "Malat1" "Fth1" "Camk2b"
## [22] "Cplx2" "Adcy1" "Meg3" "Ndrg4" "Camk2n1" "Ttyh1" "Olfm1"
## [29] "Nsf" "Plppr4"
## [1] ""
## [1] "Tmsb10" "Tuba1a" "Tubb5" "Ftl1" "Rplp0" "Ptma"
## [7] "Igfbpl1" "Rps5" "Tubb2b" "Marcksl1" "Rps9" "Tmsb4x"
## [13] "Ybx1" "Rps4x" "Eef2" "Tubb3" "Rpl10" "Nnat"
## [19] "Stmn2" "Hn1" "Ubb" "Sox11" "Prdx2" "Junb"
## [25] "Rps28" "Sox4" "Calb2" "Xist" "Gm1673" "Fxyd6"
## [1] ""
## [1] ""
## [1] "PC4"
## [1] "Tubb2b" "Tuba1a" "Zbtb20" "Tmsb10"
## [5] "Nnat" "Igfbpl1" "Tubb3" "Stmn2"
## [9] "Xist" "Marcksl1" "Tubb5" "Rtn1"
## [13] "Fabp5" "Gm1673" "Slc1a2" "Mllt11"
## [17] "Hn1" "Mt3" "Cnih2" "Sox11"
## [21] "Aldoc" "Cpe" "6330403K07Rik" "Fxyd6"
## [25] "Prdx2" "Mmd2" "Calb2" "Pantr1"
## [29] "Clu" "Cspg5"
## [1] ""
## [1] "Sparc" "Bsg" "Sepp1" "Tmsb4x" "Junb" "Arpc1b" "Hexb"
## [8] "Itm2b" "C1qa" "Fth1" "Lgmn" "C1qb" "C1qc" "Ctss"
## [15] "Ctsd" "Cldn5" "Cst3" "Jun" "Gatm" "Cd81" "P2ry12"
## [22] "Cyba" "Malat1" "B2m" "Tyrobp" "Klf2" "Trf" "Itm2a"
## [29] "Igfbp7" "Rgs10"
## [1] ""
## [1] ""
## [1] "PC5"
## [1] "Cst3" "C1qa" "Hexb" "C1qb" "C1qc" "Ctss" "Ctsd"
## [8] "P2ry12" "Lgmn" "Tyrobp" "Rgs10" "Fcrls" "Tmem119" "Csf1r"
## [15] "Fcer1g" "Selplg" "Aif1" "Junb" "Apoe" "Trem2" "Sirpa"
## [22] "Rps4x" "Ctsz" "Rps29" "Mt3" "Ctsb" "Rpl10" "Slc1a2"
## [29] "Rplp0" "Rps28"
## [1] ""
## [1] "Plp1" "Bsg" "Mbp" "Cnp" "Sept4" "Cldn5" "Cryab"
## [8] "Car2" "Gng11" "Tsc22d1" "Ptn" "Cldn11" "Mal" "Itm2a"
## [15] "Mobp" "Igfbp7" "Ptgds" "S100a16" "Mog" "Qdpr" "Aplp1"
## [22] "Pllp" "Malat1" "Mag" "Ly6c1" "Ptma" "Opalin" "Fth1"
## [29] "Olig1" "Car4"
## [1] ""
## [1] ""
qplot(x = seq(1:PCA),y = seuset@dr$pca@sdev,
xlab = "PC", ylab = "Eigenvalue")

PC25까지 사용하여 TSNE 및 Clustering
PCA_use = 25
seuset <- RunTSNE(seuset, dims.use = 1:PCA_use, do.fast = T, seed.use = 123456, perplexity=100)
seuset <- FindClusters(seuset, reduction.type="pca", dims.use = 1:PCA_use, save.SNN = TRUE, force.recalc = TRUE)
# save(seuset, file= "seuset.RData")
load("seuset.RData")
clusters_seurat = seuset@ident
# save(clusters_seurat, file= "clusters_seurat.RData")
load("clusters_seurat.RData")
Seurat 패키지 function으로 plot하기
PCAPlot(seuset)

TSNEPlot(seuset, do.label = TRUE)

ggplot2 사용하여 plot하기
# publication celltype labels
qplot(x=seuset@dr$tsne@cell.embeddings[, "tSNE_1"],
y=seuset@dr$tsne@cell.embeddings[, "tSNE_2"],
xlab ="component1", ylab = "component2",
colour = scesetFiltered$celltype)

각 cluster별 Markers 찾기
# markers <- FindAllMarkers(
# object = seuset,
# test.use = "wilcox",
# only.pos = TRUE,
# min.pct = 0.25,
# logfc.threshold = 0.25
# )
# save(markers, file="markers.RData")
load("markers.RData")
상위 5 markers heatmap 그리기
suppressPackageStartupMessages(library(dplyr))
top5 <- subset(markers, cluster %in% c(0,1,2,3,4,5,8,10)) %>% group_by(cluster) %>% do(head(.,5))
DoHeatmap(
object = seuset,
cells.use = names(seuset@ident)[seuset@ident %in% c(0,1,2,3,4,5,8,10)],
genes.use = top5$gene,
slim.col.label = TRUE,
remove.key = TRUE
)

Seurat 기능
RidgePlot(seuset, features.plot = c("Calb1"))
## Picking joint bandwidth of 0.0474

features.plot = c("Cdk1", "Ascl1", "Tfap2c", "Eomes",
"Igfbpl1","Calb2","Plk5",
"Gfap","Hes5", "Sox2", "Emx2",
"Foxg1", "Egfr", "Prom1",
"Lpar1","Nes", "Neurog2","Top2a",
"Mcm2", "Dcx", "Neurod1", "Calb1"
)
DotPlot(object = seuset, genes.plot = features.plot, plot.legend = TRUE, x.lab.rot =T)

DotPlot(object = seuset, genes.plot = unique(top5$gene), plot.legend = TRUE, x.lab.rot =T)

FeaturePlot(object = seuset, features.plot = "Igfbpl1", no.legend = FALSE,
do.hover = TRUE, data.hover = c("ident", "PC1", "nGene")
, dark.theme = TRUE)